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We study the dynamics of the g-state random bond Potts ferromagnet on the square lattice at 
its critical point by Monte Carlo simulations with single spin-flip dynamics. We concentrate on 
q = 3 and q = 24 and find, in both cases, conventional, rather than activated, dynamics. We also 
look at the distribution of relaxation times among different samples, finding different results for 
the two q values. For q — 3 the relative variance of the relaxation time r at the critical point is 
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has a finite relative variance. We speculate that this difference occurs because the transition of the 
corresponding pure system is second order for q = 3 but first order for q — 24. 
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I. INTRODUCTION 
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The random g-state Potts modclO in two-dimensions 
has recently been the subject of extensive study. One rea- 
son for this interest is that disorder even changes the or- 
der of the transition for q > 4, whejje.it is first order in the 
pure case but must be continuousotju in two-dimensions 
for all-g. Another reason interest is the prediction Jw 
Ludwigq, later confirmed by numerical simulations&QEl, 
of multi- fractal exponents at the critical point. Although 
there has beerL-a-lat.Qf. numerical work on the static crit- 
ical exponentsM'ME#jMy, there has been less on 
the dynamics, though there have been smpe-cecent stud- 
ies based on the "short-time dynamics" E^rEjO approach. 
Here we perform a careful analysis of the dynamics of the 
Potts model in two dimensions by Monte Carlo simula- 
tions. Our motivation is two-fold: 

Firstly, CardyElhas raised the possibility that, at least 
for certain types of dynamics, one might have activated 
dynamical scaling, in which the log of the characteristic 
time r (i.e. the barrier height) varies with a power of an 
appropriate length scale I i.e. lrtT ~ ft where -0 is the 
"barrier" exponent. This is in contrast to conventional 
dynamical scaling in which r itself varies with a power 
of I. i.e. t ~ l z , where z is the dynamical exponent. Ac- 
tivated dynamics has been proposed by Fisheita for the 
random field Ising model where the fixed point (of the 
renormalization group) is at T = 0. First order transi- 
tions can also be described by a T — fixed point and 
also have activated dynamics. In that case the barrier ex- 
ponent tp is equal to d — 1 for a discrete broken symmetry, 
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where d is the dimension, so tp = 1 here. 

We study two values of q: q = 3 where the pure system 
has a second order transition and so activated dynamics 
seems unlikely for the random case, and q — 24 where 
the pure system has a first order transition so there is 
a greater possibility that the random problem will have 
activated dynamics. The size of the first order jump in- 
creases as q gets large, and we find that we need q > 16 
to clearly see activated dynamics for the pure system for 
the range of sizes studied. Hence, for the random case, 
we take a rather large number of Potts states, q = 24, 
since activated dynamics is easily seen in this case for the 
pure system, so we feel we have a good chance of seeing 
it in the random case too, if it occurs. In fact, we find 
conventional dynamical scaling for both values of q, but, 
as we shall discuss, have some concern as whether this is 
the true asymptotic behavior for q = 24. 

Secondly, there has recently been considerable discus- 
sion about lack of self averaging for static quantitie 
random systems at the critical point. It is foundEJt 
that if disorder is irrelevant, then the relative variance, 
Rx, of a static quantity X tends to zero in the thermo- 
dynamic limit like L a ' v , where L is the lattice size, and 
a (< 0) and v are the specific heat and correlation length 
exponents, respectively. This is called weak self averag- 
ing. However, if disorder is relevant then Rx tends to a 
universal constant, indicating lack of self-averaging. In 
this paper we consider the question of self averaging for 
dynamics, which has not been discussed before, to our 
knowledge. Disorder is relevant for both values of q that 
we study, and yet we find rather different distributions of 
the relaxation time r in the two cases. For q = 3 we find 
a finite relative variance for r but, perhaps surprisingly, 
for q = 24 this appears to diverge in the thermodynamic 
limit and it is lnr which has a finite relative variance. 

Sec. II describes the model and the quantities we cal- 
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culate. Results for the pure models are given in Sec. 
while results for the random system are given in Sec. 
Our conclusions are summarized in Sec. |v|. 



by 



II. THE MODEL 



The Hamiltonian of the g-state Potts model is given 



(1) 



where each site i on an N — L x L square lattice is in one 
of q-states, characterized by an integer — 1, 2, • • • , q. 
The couplings, i£y, are positive, and include the factor 
of (3 = 1/fceT. They are independent random variables, 
drawn from a probability distribution, P(K). 

For the pure case, the systefl^-js at criticality if the 
"dual coupling" K*, defined by 
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K* 
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is equal to K, i.e. K* = K = K c , where e Kc = 1 + ^fq. 
For the random case it is possible to choose a self-dual 
distribution of the couplings to ensure that the system is 
at the critical point. pWe take the distribution suggested 
by Olson and Youngta (OY) since this is both very broad 
(so the system is far from the pure fixed point), and also 
does not have a substantial weight near K = (so the 
system is far from the percolation fixed point). In terms 
the OY distribution isa 



of x = e K 



Px{x) 



7T (1 — x) 2 + qx 2 ' 



(3) 



To generate random numbers with probability Px (x) one 
takes x to be 
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where r is a random number with a uniform distribution 
between and 1. We expect that the same results would 
be obtained asymptotically for any reasonable distribu- 
tion. However, to verify this would require a very large 
computational effort, since other distributions are likely 
to have larger corrections to scaling than the OY distri- 
bution, and so would need simulations on larger lattices 
than we have been able to study here. We have some 
preliminary results for a binary distribution which are 
similar to those presented here for the OY distribution, 
but a much larger computational effort would be needed 
to verify convincingly that the two distributions are in- 
deed in the same universality class. 

We focus on the time dependent magnetization 
squared, defined by 



m{tf 
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where p n (t) is the fraction of sites in state n at time 
t. The initial configuration is completely random and so 
m(t) 2 increases, eventually reaching its equilibrium value 
for t — > oo. Note that m{t) 2 is invariant under global 
symmetry transformations which permute the states n. 
Particularly useful is the average value of m(t) 2 normal- 
ized by its equilibrium value m(oo) 2 , i.e. 



m{tf 



[<m(co)2)] £ 
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where [• • -] av denotes an average over disorder, and (■ • •) 
denotes a thermal average. Clearly m(0) 2 ~ (for finitc- 
N it is 0(1/N)), and lim^oo m(t) 2 = 1. 

In addition to results for the average decay of the mag- 
netization, we also discuss the distribution of relaxation 
times of the total magnetization. 

Since we want to use a realistic form for the dynamics, 
we use standard Monte Carlo methods-tather than one of 
the more efficient cluster algorithmsE3Ej that reduce crit- 
ical slowing down. We expect that any "local" dynamics 
would give the same results, but that cluster algorithms, 
in which the average size of the cluster of flipped spins 
diverges at the critical point, would yield faster dynam- 
ics, as indeed they are designed to do. We therefore used 
the Wolfic£l cluster algorithm to obtain more accurate es- 
timates for some of the equilibrium values of m 2 , needed 
in the denominator of Eq. (H). 



III. RESULTS FOR THE PURE SYSTEM 

The transition is second order for q < 4, so we expect 
conventional dynamical scaling in this region. To see 
this we show in Fig. [I] data for fh(t) 2 for q — 3 against 
ln(t/Ti), where tl is chosen for each L in order to collapse 
the data. A log-log plot of Tl against L is shown in 
the inset to Fig. 0, which demonstrates nice power-law 
scaling. A fit, omitting the L = 4 data point, gives a 
dynamical exponent z = 2.18±0.04. An activated scaling 
plot of In tl against does not work unless ip is chosen 
to be extremely small, in which case this is equivalent to 
a conventional scaling plot. 

By contrast, for q > 4 the transition is first order which 
leads to activated scaling with ip = 1, as discussed in 
Sec. I. The correlation length at the critical point has 
been computed in Ref. |2| and we show a table of results, 
computed from the expressions in this reference, for cer- 
tain values of q in Table [jj In order to see the first order 
nature of the transition in numerical simulations the sizes 
studied must be larger than the critical correlation length 
which suggests that we should take q > 16. 

Data for q = 24 for m(t) 2 is shown in Fig. ||. The 
data collapse is good except for very short times, which 
is presumably not in the scaling regime. Interestingly, 
the data for q — 3 in Fig. [|, does scale even at very short 
times. 

Activated and conventional scaling plots for q = 24 are 
shown in Figs. [}]and [^respectively. The activated scaling 
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FIG. 1: Results for fh(t) 2 , defined in Eqs. @ and (|), for the 
pure system for q = 3. The horizontal axis is t/rx, where tl 
is determined for each lattice size by requiring that the data 
collapse. The inset shows the resulting values for tl- Since all 
the tl can be multiplied by the same constant with no effect 
on the quality of the data collapse, we arbitrarily set T4 = 1. 
The best fit, omitting the L — 4 point, gives z = 2.18 ± 0.04. 



q 


^crit 


4 


OO 


5 


2512 


6 


158.9 


8 


23.88 


12 


6.548 


16 


3.746 


24 


2.155 


32 


1.608 



TABLE I: The values of the correlation length at the critical 
point of the pure system for certain values of q. These results 
are obtained from expressions in Ref. 

fit works well, but with a barrier exponent tf> ~ 0.7±0.1, 
rather than the expected value of 1. The inset shows 
the data plotted with -0 = 1. The fit has a probability 
(Q-factor) of 6.1 x 1CP 12 , which is very low, whereas the 
fit with tp — 0.7 has a Q-factor of 0.92, which is good. 
Presumably the value ip = 0.7 is only an effective expo- 
nent which fits the data for the range of sizes studied, 
since one expects to see tp — 1 for sufficiently large sizes. 
As shown in Table |[ the correlation length at criticality 
decreases with increasing q so one expects to be closer to 



FIG. 2: Results for rh{t) 2 for the pure system for q = 24. The 
horizontal axis t/rL, where tl is determined for each lattice 
size by requiring that the data collapses as well as ^possible. 
The resulting values for the tl are plotted in Figs. |^ and 
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FIG. 3: A plot of tl (on a log scale) against 1/ for the 
pure system for q — 24. The data fits a straight line well, 
indicating activated scaling, though fitting gives t[> = 0.7±0.1, 
rather than the asymptotic value of 1. The inset shows a plot 
assuming ij) = 1. Distinct curvature indicates a less good fit 
as discussed in the text. 
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FIG. 4: A log-log plot of tl against L for the pure system 
for q = 24. The data shows pronounced curvature, indicating 
that conventional dynamical scaling does not work. 



FIG. 5: Results for m(t) 2 for the random system with the 
Olson- Young distribution, Eq. (|^), for q = 24. The horizon- 
tal axis t/rL, where tl is determined for each lattice size by 
requiring that the data collapses as well as possible. Results 
for the tl are shown in the inset and Fig. ^. The inset shows 
a log-log plot of tl against L. The straight line fit works well 
and gives a slope of z = 3.76 ± 0.04. 



the asymptotic value of ip at large q for the modest range 
of sizes that we can simulate. Our results are consis- 
tent with this since we find effective values of if) equal to 
0.55,0.7 and 0.9 for q = 16,24 and 32, respectively. For 
q — 8, activated scaling only worked for if> around 0.1 or 
less, which is sufficiently small that it is not significantly 
different from conventional scaling. Thus, although one 
expects activated scaling asymptotically for q — 8, one 
is far from this regime for the range of sizes that we can 
study. This is not surprising since the the correlation 
length at criticality, shown in Table |, is about 24 lattice 
spacings. Recently, Ozoguz et al.E3 have studied the dy- 
namics of the pure Potts paodel at criticality for q = 6 
and 7. By using the WolfO algorithm, they are able to 
study larger sizes than us, and they also incorporate cor- 
rections to finite size scaling. In this way, they find that 
their results are consistent with activated scaling with 
ij) = l. 



The main conclusion from this section is that activated 
dynamics is only seen easily in the pure model for q > 16. 
To determine if the random model has activated dynam- 
ics, we expect that one should choose q to be larger than 
the value where it can be seen for the pure system, i.e. 
we need q > 16. 



IV. RESULTS FOR THE RANDOM SYSTEM 

A. Dynamics of the averaged order parameter 

Data for m(t) 2 for the random system with the Olson- 
Young distribution, Eq. (0) , for q — 24 is shown in Fig. [|. 
The relaxation times, tj, have been determined by re- 
quiring that the data for large times collapses, (with n 
arbitrarily set to unity), and are shown on a log- log plot 
in the inset to Fig [j| The fit is good and the slope gives 
the dynamical exponent, 



3.76 ±0.04 



(9 = 24). 



(7) 



We have also attempted to scale the data using acti- 
vated dynamical scaling. However, we find only satisfac- 
tory fits are for if) very small, see Fig. [| which are not 
significantly different from the power law fit in the inset 
to Fig. H The inset to Fig. ^| shows a plot with a larger 
value of if>, 0.7. The large curvature indicates that this 
does not work. 

It appears, then, that the random Potts model has con- 
ventional dynamical scaling. However, when we discuss 
distributions of relaxations times in Sec. IV B we shall 



see that the situation is rather more complicated, and 
the true asymptotic behavior is perhaps not clear. 
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FIG. 6: An activated scaling plot of tl for the random 
system with the Olson- Young distribution, for q — 24. The 
only fits which work well have very small ip, if) = 0.1 is shown 
here, which, for the range of sizes studied, is not significantly 
different from the power law fit shown in the inset to Fig. Bl 
The inset shows a plot wiihjb = 0.7, which is the best value 
for the pure system, see Fig. H, but which clearly doesn't work 
for the random case because of the pronounced curvature. 



We have also investigated the case of q = 3, which is 
unlikely to have activated dynamical scaling, since the 
pure system has a continuous transition, and indeed we 
find that conventional scaling works well, see Fig. |. The 
fit gives 



z = 3.24 ±0.03 



(9 = 3). 



(8) 



From Eqs. (g) and (|) it appears that z increases with 
increasing q and so our results are compatible with the 
value 3.41 ± 0.06 obtained by Pan et al.l— 1 for q = 8. 



B. Distribution of relaxation times 

Since disorders-relevant for all q greater than two, 
it is expectedndO'EHl that static quantities are not self- 
averaging at the critical point. Presumably dynami- 
cal quantities like the relaxation time are also not self- 
averaging but this does not seem to have been much dis- 
cussed up to now. We have therefore also studied the 
distribution of relaxation times for the q = 3 and q = 24 
Potts models. For each realization of the disorder (i.e. 
sample), we generate the time dependent squared mag- 
netization, normalized by the equilibrium value for that 



FIG. 7: A log-log plot of tl against L for the random system 
with the Olson- Young distribution, for q — 3. The straight 
line fit works well and gives a slope of z — 3.24 it 0.03. 

sample (which we get with a good precision thanks to the 
Wolff algorithm), i.e. 



m(tf 



(m(tf) ^ 
(m(oo) 2 ) 
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The thermal average is obtained by repeating the simula- 
tion 10 4 times with the same random bonds but starting 
from different initially disordered spin configurations. 

We find that the shape of fh(t) 2 is strongly sample- 
dependent. We have tried various definitions to deter- 
mine a relaxation time r for a single sample: the time 
necessary for fh(t) 2 to reach a fixed value (0.3, 0.4, 0.5, 
0.6, 0.7, 0.8 and 0.9) and the "integrated time" defined 
by / °°(1 — fh 2 (t))dt. We checked that all these times, 
yield the same scaling behavior. The best choice, which 
minimizes the error bars, is the time to reach 0.60, which 
will be used from now on. 

We have investigated the distribution of relaxation 
times for q — 3 and 24. However, since we need to repeat 
the dynamical evolution many times (10 4 in practice), 
and since we need to repeat this for many samples, the 
range of sizes that we could study is rather restricted, 
L < 12. The number of samples used for the different 
sizes and q values is given in table Table @. 

Fig. |^ shows data for the cumulative distribution of r 
for q = 3 for different sizes with a logarithmic horizontal 
scale. One sees that the distributions for different sizes 
are shifted horizontally but otherwise look very similar. 
This is confirmed by Fig. in which the values for r are 
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FIG. 8: Cumulative distribution functions of r for the Olson- 
Young distribution of random bonds, q = 3. 



FIG. 10: Cumulative distribution functions of In t/ [In r] av for 
the Olson- Young distribution of random bonds, q = 3. This 
rescaling doesn't work since the distribution get narrower for 
larger L. 
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FIG. 9: Cumulative distribution functions of r/[r] av for the 
Olson- Young distribution of random bonds, q = 3. This 
rescaling works 



FIG. 11: Cumulative distribution functions of r for the 
Olson- Young distribution of random bonds, q — 24. 
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FIG. 12: Cumulative distribution functions of r/[r] av for the 
Olson- Young distribution of random bonds, q — 24. This 
rescaling doesn't work. 
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FIG. 13: Cumulative distribution functions of In r/[ln r] av for 
the Olson- Young distribution of random bonds, q — 24. This 
rescaling works, except for the shortest times. The inset shows 
the resulting [lnr] av plotted against InL, which gives a slope 
of zo = 3.60. 
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L 


samples 


3 


4 


2241 


3 


6 


1600 
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1537 
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12 


105 


24 


4 


2739 


24 
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2223 


24 
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868 


24 


12 


1012 



TABLE II: Parameters of the simulations used to study the 
distribution of relaxation times for the OY distribution for 
different values of q and lattice size L. 



scaled by the average relaxation time for each size. The 
data collapses moderately well, indicating that r/[r] av 
has a (non-trivial) distribution which is independent of 
size, and which is presumably also universal. Hence all 
powers of moments of the form [r™]av , (and the corre- 
sponding cumulants) all vary with L in the same way as 
the mean, i.e. as L z . In particular, the relative variance, 



R T = 



\t] 2 

I ' lav 



(10) 



is a constant for L — > oo for q = 3. This is analogous to 
the lack of self-averaging in static quantiti es .. a1. | th e crit- 
ical point which has been discussed beforc£3Elln3. Note 
that the distribution of lnr/flnrjav is sharp at large L, 
i.e. lnr/[lnr] av is self-averaging. We show the sharpen- 
ing of the distribution of lnT/[lnr] av in Fig. |l0| . 

The average relaxation time, [r] av , defined in this sub- 
section is not precisely the same as the relaxation time tl 
defined in terms of the dyna mics of the square of the av- 
erage magnetization in Sec. tVA, However, one expects 
that they should be proportional to each other, and in- 
deed this the case since a fit of ln[r] av against InL gives 
a slope of z = 3.39 ± 0.10, essentially in agreement with 
Eq. (§)_ 

Fig. |ll| shows data for the cumulative distribution of r 
for q = 24. In contrast to the data for q = 3 in Fig. [| the 
curves for larger sizes are not only shifted to the right but 
also become broader. Hence the data does not collapse 
in a plot of r/[r] av as shown in Fig. fl2"[ 

However, if wc consider the distribution of the loga- 
rithm of t then a plot of lnr/[lnr] av does scale, as shown 
in Fig. [k| Hence for q = 24, we find that R T diverges 
for L — > oo, whereas R\ nT is finite, where i?i nT is defined 
in a similar way to R T in Eq. ©. Hence for q = 24, but 
not q — 3, one can say that is is the barriers which have 
a finite relative variance. Since Fig. [l3| shows that lnr is 
the appropriate scaling variable for q = 24 we need to see 
how the average of this quantity varies with size. Hence, 
in the inset to Fig. [l^, we show a plot of [lnr] av against 
InL, which works quite well with a slope 3.60 ± 0.09. 

If we define, possibly different, exponents z n by 
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[T™]iv n ~ L Zn then the slope in the inset to Fig. is 
zq. Since [lnr] av ~ InL we have conventional, rather 
than activated, dynamics because in the latter case one 
woul d hav e [lnr] av ~ U* . This is also what we found in 
Sec. IV A . However, since we now see that the scaling 
variable is is lnr rather than r, we need to discuss fur- 
ther the behavior of averages of powers of r. We shall 
see tha t this additional analysis implies that the results 
in Sec. IV A for q = 24 may not describe the asymptotic 
behavior for large L. 

Because lnr is the scaling variable, we can write 
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In : 



dlii', 



[In r] av y [lnr], 
i nx ^ aT ^f(x) dx 

i nx< ~ a+ZolnLS >f(x)dx, 



(11) 



where f(x) is the scaling function for x = lnr/[lnr] av , 
and the last equality has used the fit shown in the inset 
to Fig. E| 

From Eq. ( pi] ) we see that the behavior of [r"] av for 
large L depends on the form of the scaling function f{x) 
for large x. If there is a very sharp cutoff at x* say, 
then the integral is dominated by values in the vicinity 

Of the Cutoff and SO [r"] av - e nx*(a+z lnL) ^ L nx*z ^ 

which gives z n = x*zq, independent of n, for n > 0. 
This is conventional dynamical scaling, except that the 
value of zq is different from that of z n with n > 0. How- 
ever, more generally, [r n ]av would not vary as L z ™, and 
hence would not correspond to conventional dynamical 
scaling. For example, a Gaussian form for f(x) would 
give [r n ] av ~ exp(const. (lni) 2 ), which increases with 
L faster than than any power but is slower than expo- 
nential. This behavior is in between conventional and 
activated dynamics. 

We have tried to estimate the behavior of [r"] av for 
large L by noting that 



where 



: g(-m[lnr] av ), 



g(k)= / e^f{x)dx, 



(12) 



(13) 



is the Fourier transform of the scaling function fix). 
Since f(x) is normalized and has mean unity by defi- 
nition, we can write 



g(k) = exp(i/c + u(k)), 



where 



*(*) = E 



(x n ) c (ik) n 



n=2 



III 



(14) 



(15) 



in which (x n ) c denotes the cumulant average of x n . 
Hence, if we can determine the analytic form of u(k) for 
large k, we can analytically continue to the imaginary 
axis and thereby obtain [r"] av 
Analyzing the data in Fig. 
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we find that the first 
few cumulants of f(x) are roughly of the form (—0.13)™. 
If we use this form then 

u(k) = exp(-0.13ifc) - 1 + 0.13ift, (16) 

which, from Eqs. ( |l2| ) and (|l4|), gives for L — ► oo 

[r n ]av ~ exp(1.13nzolnL) = L 1A3nz °, (17) 

where we have used the fit in the inset to Fig. [l3]. This 
corresponds to the "sharp cut-off" mentioned above, i.e. 
one has dynamical scaling (except that zo is different 
from z n for n > 0). However, the data is not good enough 
to be able to make this analytic continuation with any 
confidence. 

As a result of this analysis of the distribution of lnr, 
we infer that the value of z for q = 24 found above from 
b~]av (and from tl in Sec. IV A , see the inset to Fig. ||) 
may only represent an effective exponent, valid for fairly 
small sizes. 



V. CONCLUSIONS 

We have studied the dynamics of the g-state random 
bond Potts ferromagnet at the critical point in two- 
dimensions for q = 3 and q — 24. In both cases we find 
conventional dynamical scaling, and our results for the 
dynamical exponent are compatible with those of Ref . |l6| . 
However, for q — 3, we find that the reduced variance of 
the relaxation time r tends to a finite value for large L, 
while for q — 24 it is the reduced variance of In r which 
tends to a finite value. Since lnr is proportional to a 
barrier height, our results for q = 24 imply that there 
is scaling for the distribution of barrier heights but the 
characteristic barrier height scales with the log of the 
system size. 

We do not have an intuitive explanation for this dif- 
ference in behavior between the two cases, but it may be 
related to the transition of the pure system being first 
order for q — 24 and second order for q = 3. It is cur- 
rently unclear whether the scaling of the distribution of 
barrier heights for q — 24 implies that dynamical scaling 
is actually activated for L — > oo, in which case our finite 
value of z would just be an effective exponent valid for 
rather small sizes. 
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